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In this communication, the convergence of the 1/t and Wang - Landau algorithms in the calcula- 
tion of multidimensional numerical integrals is analyzed. Both simulation methods are applied to a 
wide variety of integrals without restrictions in one, two and higher dimensions. The errors between 
the exact and the calculated values of the integral are obtained and the efficiency and accuracy of the 
methods are determined by their dynamical behavior. The comparison between both methods and 
the simple sampling Monte Carlo method is also reported. It is observed that the time dependence 
of the errors calculated with 1/t algorithm goes as N -1 ^ 2 (with N the MC trials) in quantitative 
agreement with the simple sampling Monte Carlo method. It is also showed that the error for the 
Wang - Landau algorithm saturates in time evidencing the non-convergence of the methods. The 
sources for the error are also determined. 
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I. INTRODUCTION 

It is well-known that Wang-Landau (WL) algorithm is 
one of the most refreshing variations of the Monte Carlo 
simulation methods introduced in the last time Its 
effectiveness is based on the simplicity and versatility of 
the algorithm to calculate the density of state g(E) with 
high accuracy (here g(E) represents the number of all 
possible states or configurations for an energy level E 
of a given physical system). In fact, on visiting states 
with energy E the running estimate g e (E) is multiplied 
by the refinement parameter / > 1, which forces the 
system to visit less explored energy regions through the 
bias acceptance probability of min[l, —f^y ] for a move 
Ei — > Ef (with i and / the initial and final states) and 
enables a fast performance compared to other flat energy 
histogram methods |2|. 

The algorithm has been successfully used in many 
problems of statistical physics, bi ophy sics and others 
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EliEalialSl^lilli^The method is based on inde- 
pendent random walks which are performed over adjacent 
overlapped energy regions, providing the density of states 
[J, EH- In that way, thermodynamic observables, includ- 
ing the free energy over a wide range of temperature, can 
be calculated with one single simulation. 

There have been several papers in recent years dealing 
with improvements and sophisticated implementation of 
the WL iterative process in discrete and continuous sys- 
tems |, i El m EE EE EE Ei, EE EE EH, E3, HI. 

However one of the most controversial point in the appli- 
cation of the WL and others variations of the algorithm 
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is the saturation of the error between the calculated and 
the real g(E). In fact, its approaches to a constant as a 
function of the MC time, for enough long time. 

This problem was firstly evidenced by Q. Yang and J. 
J. de Pablo in reference 12611 . On the other hand, other 
authors [H, El EE, El, E3| have studied the accuracy, 
efficiency and convergence of WL algorithm. Some of 
them [H, EE E3 have demonstrated the convergence of 
WL algorithm by different arguments. Particularly, Zhou 
and Batt L 28y have present a mathematical analysis of 
the WL algorithm. They give a proof of the convergence 
and the sources of errors for the WL algorithm and the 
strategies for improvement. 

The saturation of the error and the convergence of the 
algorithm are certainly two contradictory results. More- 
over, is a crucial breakpoint because the saturation of 
the error means the non convergence of the al gori thm. 
This polemic point has been treated in reference |38l. E^j , 
where an analytical demonstration of the non conver- 
gence in the original version of the WL algorithm has 
been presented [39']. Alternatively, the authors have de- 
duced analytically the way to avoid the saturation of 
the error given an adequate form to the refinement pa- 
rameter. In fact, those methods in which the refine- 
ment parameter vary lower faster than 1/t (with t the 
Monte Carlo time) determine that the calculated density 
of states reaches a constant value for long times, therefore 
the error saturates. To overcome this limitation, they in- 
troduced a modified algorithm in which the refinement 
parameter is scaled down as 1/t instead of exponentially 
[38l |. This new algorithm allows the calculation of the 
density of states faster and more accurately than with the 
original WL algorithm due to the fact that the calculated 
density of states function approaches asymptotically the 
exact value. The 1/t algorithm has been successfully ap- 
plied to several statistical system [H, EH including the 
protein folding [40T | . 

The non convergence of the original WL algorithm 
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and other previous version, including N-fold way method 
[4ll [42j], seemed very difficult to believe. However, is 
interesting to emphasize that Landau and co-workers in 
Ref. [HI suggest that ln(ffi na i) cannot be chosen arbi- 
trarily small or the modified ln[g{E)\ will not differ from 
the unmodified one to within the number of digits in the 
double precision numbers used in the simulation. If this 
happens, the algorithm no longer converges to the true 
value, and the program may run forever. If ln(ffi na i) is 
within the double precision range is too small, the calcu- 
lation might take excessively long to finish. 

Although, the saturation of the error and consequently 
the non convergence of the WL algorithm have been 
demonstrated for a discrete system [39| , namely the Ising 
Model. The mathematical arguments of the source of 
the error for the WL algorithm seem to be more general 
and can be extended to all algorithms which consider a 
refinement parameter that change, according to the flat- 
ness condition of the energy histogram, with a law that 
decreases faster than 1/t. In all these cases, a saturation 
of the error for the calculation of the density of states and 
consequently the non convergence of the methods can be 
guaranteed. 

Recently, Y. W. Li et al. [44| report a new application 
of the well-known Wang-Landau algorithm sampling to 
the simplest continuous systems, namely the numerical 
integration. The basic idea of this new application of 
the WL algorithm is to establish a parallel between the 
density of states g(E) and the distribution g(y). Here 
g(y) represents the fraction of the integration domain 
([a, b] in one-dimension) that lies within a certain interval 
[y,y + dy]. 

This idea was proposed firstly by Troter and Dellago 
[29l ]. The authors adapted a Wang-Landau sampling 
scheme to the problem of numerical integration as an 
application of their self-adaptive range Wang-Landau al- 
gorithm. In their approach, the integrand y(x) is ex- 
pressed in terms of a "Boltzmann factor" e - "^' with 
4>(x) = ln(y(x)) and fc^T = 1 and a random walk in 
this so defined energy space is performed. Simultane- 
ously, Liang [3(| was developed a generalization of the 
WL algorithm to continuous systems. This methods was 
used to the MC integration and MC optimization. The 
scheme facilitates the numerical integration in case of 
sharply peaked functions. However, both methods ex- 
hibit a severe restriction, namely y(x) > 0. 

In principle, the WL method of integration [44| 
presents various advantages on the conventional MC in- 
tegration scheme. Specifically, these are: (i) it provides 
a procedure for the numerical integration of sharply- 
peaked or ill-behaved integrand which is difficult to be 
dealt with conventional MC methods, (ii) It can be used 
for integrands with negative values. The correspondence 
between the density of states g(E) (physical system) and 
g(y) (integration) does not require a functional form as is 
proposed in reference [1^, [3(| • (hi) It is not necessary to 
known the boundaries of the integrand such as the global 
minimum and maximum of the function y(x) within the 



integration domain, iv) It seems that the flatness crite- 
rion p and the bin width dy provide two adjustable pa- 
rameters which allow the control of the accuracy of the 
numerical estimate. The bin width seems to be the pre- 
dominant parameter for attaining reasonable estimates. 

Although the accuracy of the WL algorithm is worse 
than the simple sampling MC, at least for one and 
two dimension, as is shown in ref. its potential 

rather comes up for ill-behaved integrals and for higher- 
dimensional integration problems in general since the 
random walk remains one-dimensional. However, the 
problem of convergence of the algorithm is of a crucial 
importance in the accuracy of the numerical calculation 
for certain integrals. 

Moreover, the multidimensional numerical integration 
seems to be the more adequate test laboratory to prove 
the convergence of the WL and any other algorithms. 
There are two main reasons for that, the first one is asso- 
ciate to the fact that, as a difference with the statistical 
models where the real g(E) is only known in few discrete 
cases, for numerical integration real g(y) can be calcu- 
lated exactly for some one-to-one functions. On the other 
hand, for those well-behaved functions, it is always pos- 
sible to calculate the numerical integral using the simple 
sample Monte Carlo simulation. 

Beside that, for continuous systems, there are few sim- 
ulations where comparison has been made with the ex- 
act density of states (45[. The reason for that may be 
attributed to the non-availability of results of exact cal- 
culation for any non-trivial system having a continuous 
energy spectrum. 

In order to prove such arguments, an analysis of the 
convergence and saturation of the error for both, the WL 
and 1/t algorithms is presented. The study is developed 
in the framework of numerical calculation of multidimen- 
sional integrals. 

It is necessary to emphasize that our objective is to 
discriminate the source of errors in both algorithms by 
using the dynamical behavior of the error, instead of pre- 
senting a new method of multidimensional integration. 

The outline of the present paper is as follows: in Sec- 
tion 2, it is introduced the 1/t algorithm adapted to the 
numerical integration. In Section 3, several examples 
which include up to six-dimensions numerical integrals 
are introduced. The sources of the errors in both meth- 
ods are discussed. Finally, the conclusions are given in 
Section 4. 



II. 1/t ALGORITHM AND THE NUMERICAL 
INTEGRATION 

It is well-known that Monte Carlo methods are an 
efficient alternative to calculate numerical integrals in 
higher dimensions. Technically, Monte Carlo integra- 
tion is numerical quadrature using pseudorandom num- 
bers [H, S3, EE El- That is, Monte Carlo integration 
methods are algorithms for the approximate evaluation 
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of definite integrals, usually multidimensional ones. The 
usual algorithms evaluate the integrand at a regular grid. 
Monte Carlo methods, however, randomly choose the 
points at which the integrand is evaluated. 

The first versions of the method are rather limited, for 
instance, simple sampling Monte Carlo integration suf- 
fers from slow convergence requiring a large amount of 
sampling to reduce the statistical error. However, con- 
vergence is even not always assured. 

Several variations of the classical MC integration have 
been introduced in order to improve the performance of 
the method. 

For instance, the importance sampling Monte Carlo 
method [46| reduces considerably the statistical error 
by sampling points which are generated according to 
a probability distribution p(x) and the flattened ratio 
y(x)/p(x) is integrated instead of the original integrand 
y(x). The main limitation of this procedure arises be- 
cause the weighting probability has to be positive and 
normalized to unity in the integration domain [4q |. On 
the other hand, importance sampling methods may even 
converge to incorrect values if a bad weighting function 
is chosen; however, such errors are not readily detected. 
Moreover, conventional MC integration methods fail or 
are less efficient in case of sharply peaked or ill-behaved 
functions on multidimensional domains. 

Other methods have been introduced in the past to im- 
prove the importance sampling Monte Carlo integration, 
as for example the VEGAS algorithm It samples 

points from the probability distribution described by the 
function |/| , so that the points are concentrated in the 
regions that make the largest contribution to the integral. 

In this section, the 1/t algorithm is adapted to the 
numerical calculation of multidimensional integrals. The 
basic idea is as follows: To evaluate the definite inte- 
gral J y(x)dx is necessary to determine the proportion 
of integration domain that lies within a certain interval 
[y, y + dy], i.e. the measure {x \ x £ [a, b], y < y(x) < 
y + dy}. The distribution g{y) can be generated measur- 
ing this fraction. As it point out below, this quantity is 
a direct analogy to the density of states g(E) of a physi- 
cal system. Provided that the lower bound y m in and the 
upper bound y m ax of the integral are known the integral 
can then be approximated by 



/ = 



/ y(x)dx w ^2 9(y)y 



(i) 



To build the distribution g(y), the interval [y ma x — 
y mm ] is divided in L = [y max - y m in]/dy segment. The 
MC time is defined as t — N/L, where N is the number 
of Monte Carlo trials. In what follows all the quantities 
will be related to the numbers of MC trials N, in order 
to compare the algorithms. 

In practice, the relation S(y) = ln[g(y)] is generally 
used, in order to fit all possible values of g(y) into double 
precision numbers. 

The algorithm is as follows: 



i) Choose a value of Xi at random and the correspond- 
ing value of yi is calculated; then set S(y) — ln[g{y)\ = 
for all value of y, F a = 1 and fix Ffi na i or equivalently 

t final — 

ii) A value xj is also chosen at random and the system 
changes from yi to yf according to the probability given 
by 



P{y% -> yf) = min \ 1 



min 1 1, e 



g(ui) | 
g(yf) J 



(2) 



iii) Increment S(y) — > S(y) + F k . 

iv) After some fixed sweeps (i.e., 1000 MC time) check 
that all the sites "y" corresponding to the same Fk will 
be visited by the random walker at least one time, then 
refine F k = F k /2. 

v) If F k+1 <l/t = L/N then do F k+1 = F(t) =l/t = 
L/N. In what follows F(t) is updated at each MC time. 
The step iv) is not used for the rest of the experiment. 

vi) If t > t fi na i (F(t) < F fina i) then the process is 
stopped. Otherwise go to ii). 

Note that, with the exception of S(y), the 1/t algo- 
rithm does not use any auxiliary histogram in the calcu- 
lation of distribution g(y). 

At short times, the random walker must visit a given 
site y at least one time with the same F k . As soon as the 
refinement parameter takes F = 1/t functionality, it goes 
down independently of the number of times that the site 
y is visited by the random walker. 

As in the original WL procedure [l|, |44| , our algorithm 
provides only a relative distribution function g(y); how- 
ever, in order to evaluate the integral, g(y) needs to be 
normalized appropriately. In one-dimension, the normal- 
ized g n orm{y) is obtained by 



(b-a)g(y) 
ElZlaiy) 



(3) 



The lower and the upper bounds, y m in and y m ax, re- 
spectively, of the integrand y(x), as well as y values that 
cannot be reached within the integration domain, have 
to be determined beforehand in order to ensure the fea- 
sibility of the procedure. As is discussed previously in 
reference [2!| , one possible way to find the valid range in 
y-space is to carry out an initial "domain sampling run" 
with F = 1, before starting the actual iteration process. 
The 1/t algorithm can be easily generalized to higher 
dimensions as the original WL method [44| . 

The main objective of the present work is the com- 
parison of the dynamical behavior of the error using the 
1/t, Wang Landau (WL) and Simple Sampling (SS) al- 
gorithms. In what follows, the same nomenclature intro- 
duced in reference [44[ is used . Therefore the fractional 
accuracy is given by 
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a, (AT) 



Imc(N) f exact 



(4) 



where Imc(N) denotes the numerical estimate from 
Monte Carlo procedure and I eX act is the exact value of 
the integral. 



III. DISCUSSION OF THE COVERGENCE. A 
COMPARISON BETWEEN ALGORITHMS 

In this section, numerical integrals are calculated by 
using the three algorithms (1/t, WL and SS). The first 
two integrals given in ref. [44( are, 
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h D = J {x b - 4x 3 + x 2 - x) sin (4x)dx, (5) 

in one dimension, which its exact value is I\d — 
1.63564436296..., and 



I-2D = J J {x\ - X\x\ + x\x 2 + 2xi) 

x sin (4xi + 1) cos (4x2)dxidx2, (6) 

in two-dimension, which its exact values is I-id = 
-0.01797992646.... 

Then, it is evaluated the next one-dimensional integral, 

i rl 

In = 7 / V 7 ! - x 2 dx, (7) 
4 Jo 

which in the first quadrant, its exact value is 1^ = n . 
Next, it is evaluated the following multidimensional 
integrals, 

pl pi p\i=n 
InD = / / ••■ / J^J COS {iXi)dx\dx2---dXm (8) 

Jo Jo Jo i=1 
which their exact values can be easily obtain as, 

• i 1 

It is necessary to emphasize that those integrals given 
above have no particular physical or mathematical sig- 
nificance. However, their are very useful to compare the 
convergence of the three MC algorithms. 

Next, the results are shown. All the error estimates 
are obtained from 100 independent simulations. 

In the Figure 1, the fractional accuracy for the one- 
dimensional integral I\d given in eq.(5), as a function 
of the MC trials calculated by the three algorithms is 



FIG. 1: Dynamical behavior of fractional accuracy a/ for one- 
dimensional integral given in eq.(5) calculated by means of 
Wang-Landau, 1/t and simple sampling Monte Carlo integra- 
tion. The WL calculations have been made with a flatness 
criterion p — 0.9. In both cases, the WL and 1/t algorithms 
the bin width is dy — 0.005. The quantities are obtained 
averaging over 100 independent samples. 




Number of Monte Carlo trials (N) 



FIG. 2: Dynamical behavior of fractional accuracy a/ for two- 
dimensional integral given in eq.(6) calculated by means of 
Wang-Landau, 1/t and simple sampling Monte Carlo integra- 
tion. The WL calculations have been made with a flatness 
criterion p — 0.9. In both cases, the WL and 1/t algorithms 
the bin width is dy = 0.005. The quantities are obtained 
averaging over 100 independent samples. 

shown. The WL calculations have been made controlling 
the histogram every 10000 MC with a flatness criterion 
p = 0.9. In both cases, the WL and 1/t algorithms the 
bin width is dy = 0.005 and L — 3066. Clearly, the 
fractional accuracy for the WL algorithm saturates for 
N w 10 9 , while the 1/t and simple sampling calculations 
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FIG. 3: Dynamical behavior of fractional accuracy a/ for dif- 
ferent values of bin width dy by means of 1/t algorithm. The 
same integral and condition described in Figure 1 is used. 



g(y) as strategy to calculate numerical integrals. In Fig- 
ure 3, the effect of dy in the calculation of the error for the 
1 jt algorithm (for the same integral and condition given 
in Figure 1), is shown. The dependence of the error with 
the bin width will be determined by the characteristic 
of the function y(x). In fact, one can expect that for 
smaller dy ^ the saturation will occur at longer times. 
However, that is not generally valid, as is observed in the 
figure where for dy = 0.05 the error saturates before than 
for dy = 0.1. 

In Figure 4, the effect of the bin width dy on the behav- 
ior of a/ for the integral referenced in Figure 2 is shown. 
In this case the smaller values of dy lead to small errors. 
As a difference with the WL algorithm, the numbers of 
Monte Carlo trials is an input parameter for the 1/t and 
simple sampling algorithms. 

To make a deep analysis of the effect of the bin width in 
behavior of the error in the 1/t algorithm, let us consider 
the integral given in eq. (7) , which its exact value is 1^ = 
7r. For the definition of the distribution function and 
considering that f{y) is one-to-one function, it is possible 
to obtain the exact value of g(y) as, 




g ex (y) = ^l-(y + d V y - y/T 



y 



(10) 
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FIG. 4: Dynamical behavior of fractional accuracy a/ for dif- 
ferent values of bin width dy by means of 1 /t algorithm. The 
same integral and condition described in Figure 2 is used. 



of the error are in close agreement with a behavior given 
by af oc 1/y/N, as is expected. 

In Figure 2, the fractional accuracy for the two- 
dimensional integral Iid given in eq.(6), as a function 
of the MC trials calculated by the three algorithms is 
shown. In both cases, the WL and 1/t algorithms the bin 
width is dy = 0.005 and L — 1213. The flatness criterion 
for WL algorithm is the same used in Figure 1. The be- 
havior of the error is similar to the one-dimensional case 
presented in Figure 1. 

As one can observed in both cases the error for WL 
saturates, then a/ does not scale as N^ 1 ^ 2 . 

The bin width certainly introduces a systematic error 
in all the algorithms that use the distribution function 



which is valid for all dy ^ 0. Then, for a given value of dy, 
the corresponding "exact" value of the integral I 7r (dy) ex 
can be obtained. Let us define a new fractional accuracy 
related to the value of dy as 



a f (N,dy) = 



Imc(N) - I v (dy) 



(ii) 



Let us analyze the effect of the bin width in both WL 
and 1/t algorithms by calculating the integral given in 
eq.(7) and using the fractional accuracies given in eq.(4) 
and eq.(ll). 

In Figure 5 a) and b) it is shown the fractional accuracy 
defined in eq. (4) versus the MC trials for different values 
of the bin width dy for the WL and 1/t, respectively. The 
error saturates in both cases as is expected. However, in 
the next figures, it is demonstrated that the sources of the 
saturation in both cases obey different causes. In fact, to 
avoid the effect of dy on the saturation of the error, in 
Figure 6 a) and 6 b) it is shown the fractional accuracy 
defined in cq.(ll). Here for each value of dy one obtain 
the exact value of lTr(dy) ex . While for the WL results the 
error still saturates (Fig. 6 a) in 1/t calculations (Fig. 6 
b) one can observe that no saturation occurs. This is a 
clear evidence of the convergence of the 1 ft in all discrete 
systems. 

Therefore, one can conclude that the source of satura- 
tion in 1/t algorithm is exclusively due to dy. In fact, for 
a given dy the error approaches asymptotically to exact 
value of integral I^(dy) ex which is different of 7r, as one 
can see in the inset of Figure 6 a). However, for WL 
algorithm the source of the error saturation are intrin- 
sically associate to the nature of the algorithm, namely 
the decay of refinement parameter. 
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FIG. 5: Dynamical behavior of fractional accuracy a,f for one- 
dimensional integral given in eq.(7) calculated by means of a) 
Wang-Landau, and b) 1/t algorithms for different values of 
bin width dy. The WL calculations have been made with a 
flatness criterion p — 0.9. In both cases, the results are com- 
pared with the simple sampling Monte Carlo integration. All 
the quantities are obtained averaging over 100 independent 
samples. 



Finally, let us describe the behavior of both algorithms 
in the calculation of higher dimension integrals. 

In Figure 7 a) to d), it is plotted the integrals given 
in eq.(8) with n — 3,4,5,6. The numerical calculations 
have been made using the three algorithms describe in 
the paper. In all case dy = 0.05, L = 40. As one can ob- 
serve in all case the error for the WL algorithm saturates 
demonstrating that the non convergence of the method 
is independent of the dimension. On the other hand, the 
results obtained by using the 1/t algorithm is in close 
agreement with simple sampling. 

The numerical estimates of statistical error for all the 
integrals are given in Table I, using the three algorithms 
described in the paper. The number of final Monte Carlo 
trials per run is Nfi na i = 10 10 . All the WL calculations 
have been made with a flatness criterion p = 0.9. The 
values of the estimates showed in table I confirm that 1 jt 
is more accurate than WL algorithm for all the integrals 
in any dimension. 



FIG. 6: Dynamical behavior of fractional accuracy aj for dif- 
ferent values of bin width dy, for the same integral and con- 
dition explained in Figure 5. In the inset of Figure a), the 
exact value of In(dy) ex is expressed. 



IV. CONCLUSIONS 

In this work the numerical calculations of multidimen- 
sional integrals are used to analyzed the convergence of 
1/t and WL algorithms. The numerical integration is 
an excellent laboratory to prove the convergence of the 
algorithms for different reasons: i) in many cases the in- 
tegrals can be solve analytically, then it can be easily 
checked the dynamical behavior of the error; ii) for some 
one-to-one functions, one can easily obtain the exact ex- 
pression for the distribution function, g(y). This is an 
advantage over the physical systems, where only in few 
cases the exact density of states is known, particularly 
for continuous system; ii) the initial and final states are 
not correlate with their neighborhood, namely, a given 
initial state can be changed to any other final state in 
the integration domain of y(x). In this sense, if a given 
algorithm can not converge appropriately in the calcula- 
tion of a numerical integrals, it will be more difficult to 
do it in those physical systems where the initial and final 
states are strongly correlated. 

Although the flatness criterion p and the bin width dy 
provide a source of saturation of the error in the WL 
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TABLE I: Numerical estimates of integrals calculated by using the Wang-Landau, 1/t and simple sampling Monte Carlo. The 
number of final Monte Carlo trials per run is Nf ina i — 10 10 . All the WL calculations have been made with a flatness criterion 
p = 0.9. Results and error estimates are obtained from 100 independent simulations. 



Integral WL Algorithm 1/t Algorithm Simple MC exact 

— hr> 1.635580(285)" 1.635617(27)" 1.635752(23) 1.63564436296... 

Iid -0.0179671(152)" -0.0179790(51)" -0.0179841(68) -0.01797992646... 

In 3.1415799(35) 6 3.1415819(23) 6 3.1415920(15) 3.14159265358... 

h D 0.01801608(2989) c 0.01799079(46) c 0.01799666(78) 0.01799626791... 

h D -0.00339585(1798) c -0.00340505(54) c -0.00340506(55) -0.00340490511... 

ho 0.00065747(505) c 0.00065298(47) c 0.00065225(49) 0.00065300923... 

ho -0.000031919(3288) c -0.00003003(15) c -0.00003110(38) -0.00003041015... 



"The bin width is dy = 0.005. 
6 Thc bin width is dy = 0.0005. 
The bin width is dy = 0.05. 
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FIG. 7: Dynamical behavior of fractional accuracy a/ for n- 
dimensional integral (8) with a) n = 3, b) n — 4, c) n = 5, 
and d) n = 6, using Wang-Landau, 1/t and simple sampling 
Monte Carlo integration. The WL calculations have been 
made with a flatness criterion p = 0.9. In both cases, the WL 
and 1/t algorithms the bin width is dy = 0.05. The quantities 
are obtained averaging over 100 independent samples. 



algorithm. The main reason of error saturation is that 
the refinement parameter is scaled down exponentially 
instead of a power law. 

On the other hand, it is shown that the behavior of 
the error in the WL algorithm, for a single value of the 
parameters dy and p, does not follows the 1/y/N at any 
time. 

Alternatively, the convergence of 1/t algorithm is an- 
alyzed through the calculation of multidimensional inte- 
grals. 

To obtain the density of states function (physical sys- 



tems) or the distribution function (integral calculation) 
in a continuous system, it is necessary a grid discretiza- 
tion. This introduces a systematic error in the calcula- 
tion of the observables which depend on the size of the 
unitary cell of the grid, namely the bin width. 

For this reason, the error in the 1/t algorithm saturates 
as a function of dy. However, when the continuous model 
is approached by a discrete lattice and the corresponding 
value of the distribution function can be obtained exactly, 
the calculation approaches asymptotically to the exact 
value of the integral without error saturation. 

On the other hand, the behavior of the error in 1/t 
algorithm is in close agreement with the simple sampling, 
following the 1 /\/N law. Moreover, the 1/t algorithm can 
be used as a reference in the calculation of the density of 
state in physical systems due that up to now there are 
no other method which can calculate the error below the 
limiting curve 1 / >/~N. 

In summary, in the present paper it is shown a new ev- 
idence of the saturation of the error in the WL algorithm 
which implies the non-convergence of the method. 

In contrast, the dynamical behavior of the 1/t al- 
gorithm is analyzed, concluding that the algorithm is 
always convergent in all discrete system. In the contin- 
uous model, the only source of error saturation is the 
grid discretization. Therefore 1/t algorithm is a more 
efficient, accurate and easy implemented to simulate the 
distribution function (density of state) in discrete and 
continuous systems without using any histogram. 
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